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Abstract 

We study a stochastic two-species chemical reaction system with two mechanisms. One 
mechanism consists of chemical interactions which governs the overall drift of species amounts 
in the system; the other mechanism consists of resampling, branching or splitting which makes 
unbicised perturbative changes to species amounts. Our results show that in a system with a 
large but bounded capacity, certain combinations of these two types of interactions can lead to 
stochastically-induced bistability. Depending on the relative magnitudes of the rates of these 
two sets of interactions, bistability can occur in two distinct ways with different dynamical 
signatures. 
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1 Introduction 



Recent advances in measurement technology have enabled scientists to observe molecular dynam- 
ics in single cells and to study the cell-to-cell variability (Brehm-Stecher [5]). Many studies have 
shown that variability observed in genetically identical cells is due to noise that is inherent to 
biochemical reactions happening within each cell (Arkin [20], Elowitz [10]). Understanding how 
intracellular mechanisms are affected by this intrinsic noise is an important challenge for systems 
biology. Determining what role this noise plays in creating phenotypic heterogeneity has many 
practical consequences (Avery [ ]). 

An important feature in cellular dynamics is bistability, the alternation between two different stable 
states for a molecular species. This feature is present in many gene-expression systems, where a gene 
alternates between two types of states ( "on" and "off" ) regulating the production of a protein. It is 
also present in many phopsorylation switches in signaling pathways. Causes for bistable behaviour 
can be deterministic, but many bistable switching patterns are enabled by stochastic fluctuations. 
It is often assumed that it follows from the existence of two stable equilibria in the deterministic 
drift and the ability of infrequent large fluctuations to pull the system from a basin of attraction 
of one equilibirum to the other. There are also cases of chemical dynamics in which bistability is 
not possible in the deterministic model, but is possible in the stochastic model of the same chemical 
reaction system ([ ] and [ ]). Metastable behaviour is also sometimes observed (LydRobert [24]). 

In addition to noise inherent to biochemical reactions, cells also experience fluctuations in molecular 
composition due to cell division. This source of noise is significant, and also difficult to separate 
from the noise due to biochemical reactions (Paulsson [22], [2.3]). In this paper we investigate under 
what conditions a system of chemical reactions in a cell can use these two sources of noise to exhibit 
bistable or metastable behaviour in their molecular composition. 

We would like to emphasize a couple of points observed in the literature. First, the rate of switching 
between two states is important for cellular development and survival (vanOud [1]). Time-scales 
on which transitions between stable states happen varies whole orders of magnitude over different 
systems. For example, in the lysogenic state of E. coli the time-scale of switching between states is 
slow (Zong [2!)]) - once per 10® cell generations - as determined from the activity of a controlling 
protein. In the case of gene expression in S. cerevisiae the switching time-scale is fast (Kaufmann [15]) 
- once per 8.33 generations - and switching times between mother and daughter cells are correlated 
in a way that takes several generations to dissipate. Second, both the strength and the distribution 
of noise affects whether bistability will occur and what the final outcomes will be. (Samoilov [2."^i]) 
and (Bishop [4]) show that auxiliary chemical reactions can induce a dynamic switching behaviour 
in the enzymatic PdP cycle, and that final dynamics is determined by the noise of the additional 
reactions. In the bistable switch of lactose operon of E. coli (LydRobert [24]) show that both cellular 
growth rate and the molecular concentration levels influence the ability to switch. (Paulsson [23]) 
showed that the type of the cellular division mechanism also plays an important role in the form of 
the final dynamics. We interpret these observations vis-a-vis our results in the Discussion section. 

Finally, we note that bistablility in a stochastic population system is not limited to chemical dy- 
namics. In a genetic population mutation and selection may lead to alternating fixation in one of 
two genotypes. In an ecological population interactions between species can lead to dynamics where 
two competing species are switching for dominance. We note that our analysis and results apply to 
any population model described by a density dependent Markov jump process. 
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1.1 Outline of Results 



We examine qualitatively different ways in which switching between stable states is a result of a 
stochastic effect in a population modeled by density dependent Markov jump processes. In addi- 
tion to noise inherent to the reaction system, we include an intrinsically noisy splitting/resampling 
mechanism in the system. In many stochastic branching models an entity will (upon reproduction, 
division, duplication, etc.) produce offspring identical to itself. Here we model the division as un- 
biased but variable. When a cell divides its molecules are randomly allocated to its daughter cells, 
only on average replicating the parent's molecular composition. Wc will show that introducing such 
a splitting process at a sufSciently high rate can produce switching dynamics in which previously 
unattainable states become attainable. We will exploit the fact that these two sets of mechanisms 
(reactions in the system and changes due to unbiased resampling/splitting of the system) may op- 
erate on different time-scales. 

We consider the following question: which qualitatively different types of behaviour can we observe 
and under which time scaling regimes? The short answer is as follows: (1) If the resampling mecha- 
nism is "slower" than the reaction dynamics, then the system behaviour will entirely depend on the 
nonlinear dynamics of the reactions: in case the underlying deterministic system has multiple stable 
equilibria, the stochastic process will behave as a Markov chain switching between these states. (2) If 
the resampling mechanism is much "faster" than the reaction dynamics, then the system behaviour 
will not depend on details of the reaction dynamics, and will behave as a Markov chain switching 
between two extremes (zero and capacity) of the system. We define a single parameter based on 
the rates of the two mechanisms that makes the meaning of "faster" and "slower" in the statements 
above mathematically precise. 

We show that a fast but ubiased resampling mechanism may be necessary to produce bistable 
behaviour that the reaction dynamics cannot exhibit. We further show that the two cases (1) and (2) 
produce qualitatively different dynamical signatures, in terms of switching times and stable points. 
Since our analysis only depends on general features (unbiasedness and time-scale of the rate) of the 
resampling mechanism, one can also use a set of auxiliary reactions instead of resampling. There are 
other types of noisy mechanisms that one could consider, however, our goal is to stress that adding 
noise with even small changes (relative to the size of the system) can produce bistable behaviour. 
The additional noise achieves this either by: (1) introducing small perturbations to a dynamical 
system that already has the required properties for bistability; or (2) occurring so frequently that 
the details of the dynamical system are irrelevant and the system is pushed to its extreme (zero or 
capacity) amounts. 



2 Description of the Process 

2.1 Stochastic Model for Reaction Dynamics 

In the customary notation for interaction of chemical species labeled A,B,. . . 

{aiA + hB + --- a■A + 6■B + •••}^=l,...fe (1) 

denotes a system of reactions indexed by i = 1, . . . ,k in which ai,bi,. . . e Z+ molecules of types 
A,B,... respectively react and produce a -,6-,... e Z+ molecules of these types. Each reac- 
tion i has a reaction rate A^, a time and state-dependent rate of occurrences of this reaction. If 
XA{t),XB{t), . . . denote the number of molecules of type A,B,... respectively at time t > 0, then 
X {t) — {Xa {t) , Xb (t), . . .) evolves as a Markov Jump Process with jumps sizes { (a ■ — , 6 ■ — 6^ , . . . ) } 
occurring at rates Xi{XA{t),XB{t), . . . ). 
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The reaction formalism (1) can also be used to describe other systems of interacting entities under a 
well-mixed spatial assumption. For example, evolution of an SIS epidemic is expressed as: I+S — > 21 
(infection), / — ^ 5 (recovery); a two- allele Moran model with mutation from population genetics can 
be expressed as: B,B ^ A (mutation), A + B ^ 2A, A + B ^ 2B (resamphng). 

For simplicity, we consider the effect of a system of chemical reactions on essentially a single molecular 
species A. We include the effect of only one other species B which satisfies a conservation relation 
with A. This means that every reaction involving A and B is of the form aA+bB {a+(^)A+{b-'(^)B 
for some G {—a, ...,&}, and it ensures that the state space of our system is one-dimensional 
determined by {Xa {t),t > 0) . The rationale for such a conservation law could come from a cellular 
environment which is limited (by a factor such as space, or availability of nutrients or catalysts), or 
a molecular species whose type can take two different forms (e.g. a gene that has two allelic types). 

We also assume the following properties for the reaction dynamics: 

1. The amount of species X^(<) is bounded above by the system capacity N and below by 0: the 
rate of any reaction that decreases the amount of A is zero when = 0, and the rate of any 
reaction that increases the amount of A is zero when Xa = N. 

2. The drift at and N of the overall reaction dynamics is directed towards the interior: 

^i:[XA{s)\XA{t) = 0]|^^^ > 0, ^E[XA{s)\XA{t) - N]\^^^ < 0. 

3. The form of reaction rates A is governed by the law of (stochastic) Mass- Action Kinetics: a 
reaction of the form: 

aA + bB ^ a' A + b'B 

has rate X{X{t)) ^ K{XA{t))a{XB{t))b = K{XA{t))a[N - XA{t))b- Here (Z)c denotes the 
falling factorial {Z)c = Z{Z — 1) ... (Z — c -I- 1). When we renormalize Xyi(t) by its maximum 
value N ^ we will also need the 'scaled falling factorial' {z)c.n defined by: 

(z),,Ar := iV-^(7Vz), = z (^z - (^z - • • • (^z - , 0<z<l. (2) 

Note that limjv->oo(-z)c.7V = z'^ for fixed z and c. The constant k > is independent of the 
state (X^(t),XB(t) — N ~ XA{t)) but will depend on the scaling parameter N, k = k{N). 
We do not necessarily assume that k{N) has the 'standard' scaling form n{N) = k7V^^(°+''). 

4. The effect on A from any other species in the system is subsumed into the values of the rate 
constants k, and are assumed to be state-independent. 

Assumption 1. ensures that Xa{-) G {0, . . . A^} where A^ serves as the system-size parameter, 
while Assumption 2. ensures the reaction system does not get absorbed at either boundary {0, A^}. 
Assumption 3. is not essential, but with an explicit scaling of the rate k{N) in terms of N, the 
polynomial form of the rates A will make it easy to also establish the scaling of the rates X{X{t)) in 
terms of A^ under a rescaling of the species amounts Xa (we will occasionally use the notation k{N) 
for K when awareness of dependence on A^ is key). Assumption 4. is made to absorb the effect of 
the environment and other species, the changes of which we will not keep track of explicitly. 

Under these assumptions, our reaction network system can now be expressed by: 

{aA + bB — > {a + C)A+ {b - C)B} ae{o,...,N},be{o,...,N}xe{-a,...M (3) 
with reaction rates of the form X'i''{x) ~ nf" ■ {x)a{N — x)b- 
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Since the dynamics of the system depends on its ovcraU drift, it wiU be useful to distinguish a subset 
of reactions in the system that is 'balanced', in the sense that their combined expected contribution 
to the net change in the amount of A is zero, irrespective of the current amount Xa- Let I denote 
the set of all triples (a, b, () for which a reaction as written in (3) is present in the system. A subset 
of reactions forms a 'balanced group' X^^^ C I for some reactant amounts a*, b^, if it satifies: 

C:(a.,b.,C)Gl'>'»l C:(a.,f'.,C)62:'"'' 

A reaction (a, b,C,) ^ ^ that is not part of some balanced group we call 'biased', I^^^ = I — I*^^'. 

For any balanced reaction (a,5, C) G I^^^, there is necessarily a reaction (a, 6, C') G 2^^^^' with C' 
having opposite signs (though not necessarily of the same size). Hence, a reaction (a, 6,0 G I^'^^ 
cannot have non-trivial rate at the boundaries of the system: if A"'''(0) > for some C > then 

the balance condition would imply the existence of some C < for which A^; (0) > 0, which 
would violate Assumption 1 by allowing Xa to drop below upon a single further (a, b, (') reaction. 
Consequently, the boundaries and N are absorbing for the balanced subsystem of reactions, and 
for all (a, b, (') € I^^^ we must have both a > and 6 > 0. Since Assumption 2. does not allow the 
boundary {0, N} to be absorbing for the full dynamics, this further implies that there is at least one 
biased reaction {0,b,() G I*^"* with 6, C > and A^'''(0) > 0, and also at least one biased reaction 
(a, 0, C') G I*"'^ with a > 0, C' < 0, and A^;°(Ar) < 0. 

The continuous-time Markov jump process model for the reaction dynamics can be expressed in 
terms of a set of Poisson processes under a random time change. Given a collection {y^°^}(a.6,c)e2^ 
of independent unit-rate Poisson processes, the state of the system can be expressed as a solution 
to the stochastic equation (see [10] or [■]] for details): 



XA{t)^XA{0)+ Yl ^Yt{l ^tiXA{s))ds) 

ia,b.X)€l ° 

= Xa{0) + V CY^^'i f Af (X^(s))ds) + /* F{XA{s))ds 

,„i,,^^T ^ Jo 



{a,b,0£l 

where {^'^°^}(a,6,c)ei ^''S centered Poisson processes Y{\t) := Y{\t) — \t, and 
F(x)- Y CAf(x)= Y C«f(a;)a(iV-x)b. 

(a,b,C)eI (a, 6,061'''=' 



Since the capacity N of the system may be arbitrarily large, we will consider a 'standard' rescaling 
of the system (see e.g. Ethier and Kurtz [11], Ch 11.2). Let XN{t) = N-'^XA{t), then 



XN{t) = Xn{0)+ Y1 N-^Y^-'lN'^+'K^f f iXNis))a.Nil - XN{s))b.Nds] + f FN{XN{s))ds 
(C.a,6)ei \ Jo J Jo 

(4) 

where the local drift of the renormalized system is given by 

(C,a,b)el''''' 

The most important feature of the Markov jump process model is the relationship of the variance 
to the drift. Note that we can write (4) as 

X^it) = Xn{0) + M^it) + I F^{X^{s))ds 
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where the second term from (4), a weighted sum of time-changed centered Poisson processes, is a 
martingale Mj^(t) whose quadratic variation satisfies 



iC,a,b)ex ^ 

Hence, if J^t — (j{X{s),0 < s < t) denotes the natural filtration of the process, then 

^E[Xjv(s)|J-i] =E[Fiv(Xiv(t))l-Ft] = ^ iV"+^-^C <''(Xjv(f)),,iv(l - X^v W)i,,iv 

{c,£i,6)ei''''> 



E[{XNis)~E[XNis)]f\Tt' 



ds 



^E[[A/^].|J-,] 



J2 7V'^+'-\'<'(XjvW)a,iv(l-X^(t))6,iV. 
(C,a,6)Sl'=''lul'''=' 



Recall that the reaction rates = k'^^{N) also depend on the scaling parameter N. The standard 
scaling for a reaction constant is k'^'^{N) — ii'^''N^~^°'^''^ for some TV-independent constant k^''. 
However, regardless of the chosen scaling of k"'', for biased reactions I*^'^ the order of magnitude 
for each summand in the infinitesimal variance ^E[[M7v]s|-^t] \ is times smaller than the 
corresponding summand in the infinitesimal drift E,[Fff{Xff{s))\J^t] This constrains the possible 
limiting dynamics of X^- Suppose the scaling of the rates is k^'' — A^^~'^°+'')k"^, and note that then 
Fn{x) — )■ a b)GX'=i'' C'*^''2;"(l — x)'' Uniformly for x S [0, 1]. As established in [LS], in the hmit 
as iV — >• oo the drift overpowers the noise and, provided Xjv(O) x{0), the renormalized process 
(Xjv(t),t > 0) converges in distribution (in the Skorokhod topology of cadlag paths) to a solution 
{x{t), t > 0) of the ordinary differential equation 

x{t)=xiO)+ f Ciifx{s)''{l-x{s))''ds. (5) 

In fact, if the scaling of the reaction constants k'^' is not standard, but is consistent for both 
balanced and biased reactions in terms of the polynomial order of the rate function A^'', then the 
same deterministic limit is obtained under an appropriate time rescaling. 

The only way to get a stochastic limiting object for Xn is for at least one rate constant of a balanced 
reaction to have a different scaling in N. This different scaling needs to be such that the noise term 
of this reaction will be of the same order of magnitude as the overall drift from the biased reactions. 
This would require a specific separation of time-scales for balanced versus biased reactions. Although 
we do not exclude this possibility from our analysis (see definition of ea at the end of this section), 
our emphasis in this paper is on separating the time-scales in terms of contribution of an additional 
source of noise, and its ability to produce nontrivial random limiting objects for Xj^j. 



2.2 Stochastic Model for Resampling, Branching or Splitting 

We now introduce the additional mechanism in the system that describes changes to species amounts 
due to the effect of splitting, branching or resampling, which also effects the species count. For 
intracellular molecular populations, our first model of splitting was motivated by a simple double- 
then-divide principle: the cell will first double in size by replicating its constituent molecular species, 
and then allocate approximately one half of this doubled material into each daughter cell - the 
allocation mechanism is not perfect and will make random error from the original (undoubled) 
amount. For genetic populations, common models for resampling follow the Wright-Fisher or the 
Moran neutral reproduction law: each individual of the offspring population chooses at random from 
the diploid version of the current population's genes what to inherit - the resampling mechanism 
is such that an allele of one type in one generation may at random be replaced in the subsequent 
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generation by an allele of the other type. These two are both examples of a general mechanism with 
the following key properties that we assume for splitting/resampling: 



5. The splitting/resampling occurs at rate "f{x, N) that depends on: the current state Xa = a; of 
the system and the scaling parameter N; conditional on = a; it is independent of reactions. 

6. The change in the species amount Xa due to a splitting/resampling event has the distribution: 
Px,y = P[^yi(i) — y\^A{t^) = x] that have absorbing boundaries po.o = 1,Pn,n = 1, and that 
are unbiased: 

IJ.n{x) ^^ypx,y = X y X e {0,l,...,N}. 

y 

We also assume, for some of our results, that the rate 7(2;, N) and distribution {px,y} are such that: 
T. The change sizes are asymptotically uniformly bounded: 

(7! a) VA > 0, sup 7(2;, iV) ^ Px,y ^ as N ^ 00, 

y:N-^\y-x\>A 

and the change size variance uffix) — '^y{y — x)'^px,y is asymptotically given by: 
(7rb) snp\j{x,N)N-'^a%{x)-^^a^iN~^x)\^0 as iV -> 00, 

X 

for some constant 7 > and function that are independent of N, and such that x a'^{x) 
is continuous with <t^{x) > 0,Va: G (0, 1) and o-^(O) = <t^{1) = 0. 

Unbiasedness in Assumption 6. could be replaced by an 'asymptotic unbiasedness' assumption 
N~^\iin{x) — a;| — > as iV — > 00, but for the sake of simplicity we assume hn{x) — x. Absorption 
in Assumption 6. implies splitting is noiseless on the boundaries regardless of its time-scale. When 
the additional Assumption 7* holds (as we will assume for our results in Section 2.3), the splitting 
mechanism contributes diffusively to the limit of the renormalized species count X^r. However, we 
will also examine the case when the rate of the splitting mechanism is on a slower time-scale (in 
Section 3.1), as well as the case when it is on a faster time-scale (in Section 4.1). The condition 
that has boundary values (t^(0) — (t^{1) — is natural given that any splitting or resampling 
mechanism should absorb at the boundaries as indicated by po.o — Pn,n — 1- 

Example (HG): One example of a splitting mechanism would be to completely randomly reallocate 
the doubled content of a parent cell into daughter cells. If the initial content is {Xa, Xb) — (x, N—x), 
and the doubled content (2a;, 2(A^ — x)) is partitioned in a single swoop (draw without replacement) 
into two sets of N molecules (one for each daughter cell), then the content in each daughter cell has 
the hypergeometric distribution (below we keep track of an arbitrarily chosen single lineage): 

/2x\ /2N-2x\ 

p^,y = V[XA{t) - y\XA{t-) =x\^ " '^J^'" , V (2x - iV) < 2/ < 2a; A iV. 

The change in the species count is clearly unbiased iin{x) = Y^'^^ovn-2x VPx.y — with variance 
^"'"'^ N2x{2N~2x)( N-l\_x{N-x) 



4iV2 V 2N -1 2N -1 

y = QVN-2x 



Then Assumption T. will hold if 7(3;, N) = j'^N and (j'^{x) — ^x{l — x), since for (7? a) we have 

N 

i7V(2iV-l) 2N^^ N' 



sup |7(a;, N)N-^al{x) - f5^{N-^x) \ = f sup | ^^^^ - ^ ^ (l 
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and using tail bounds for the hypergeometric distribution, [ ], ^y^x+NAPx,v < e ^'^ ^ indepen- 
dently of X, and for (7!b) we have 

7(a;,iV) sup ^ p^^y < 2^^ Ne'^^^'^ ^ 0. 

y:\y-x\>NA 

Example (Bin): Another example would be to sample with replacement from the population in 
which each offspring picks its type randomly from any individual in the parent generation. If the 
initial count is Xj^ — x, then the count in the next generation has the binomial distribution 

p.,,=P[X^(t) = y|X^(i-)=-^]= 0<y<iV. 

This form of resampling is used in (the haploid version of) the Wright-Fisher model for genetic drift 
(e.g. [ ] Sec 1.2.). It is also used as the prototype of a splitting mechanism of simple 'independent 
segregation' of division of cells [ .'"]. This distribution is again unbiased, and Assumption 7*. will 
hold if j{x, N) — \')'^N for some constant 7^ > 0. Using similar arguments as above, it is then easy 
to show that both (7*a) and (7!b) will hold with (t'^{x) = \x{\ — x). 

Example (Bern): Finally, the simplest example of a splitting/resampling mechanism is to have a 
single amount error in the daughter cell (or the next generation), and to have the rate at which the 
error occurs be proportional to both the current amount = x and the amount of Xb = N — x. 
Errors from imperfect division will result in ± change with equal probability 

Px,x-1 = Px,x+1 = 1/2- 

This distribution is clearly unbiased, and Assumption 7* will hold if the rate of error occurrences is 
7(2;, N) = ^^^N'^jj{1 — j^) for some 7^ > 0, with the limiting variance a^{x) = ^x{l — x). 
This form of resampling is used in the Moran model for genetic drift (e.g. [ ] Sec 1.5.). It is also used 
in [2.')] as an example of an 'ordered segregation' splitting mechanism for cell division (self volume 
exclusion partitioning error, [ ] Supporting Information p5). In a cellular system it could also be 
described as a set of balanced reactions A + B ^ 2A, A + B ^ 2B with mass-action dynamics and 
appropriately scaled rate constants. 

We note that, from the perspective of limiting results, the differences in the specific details of the 
mechanism will not be important. The only feature of relevance will be the order of magnitude of 
the prelimiting rate j{x, N) and the form of the limiting variance a^{x). There are many other types 
of splitting, branching or resampling mechanisms, yielding a different form for the limiting variance. 
They are easy to construct in case of small changes that result in single count errors, via a range of 
birth-death probability distributions. We shall see, in both Section 3 and Section 4, how the actual 
form for the variance a'^{x) affects the qualitative behaviour of the limit of the renormalized process. 

The changes due to this additional mechanism can also be expressed in terms of a Poisson processes 
under a random time change. Let 1^ be a counting process with state-dependent rate 7(0;, iV), and 
{Z{x, s)}q<cx<n be independent random variables with probability distribution p^^. for any s > 0. A 
change due to splitting or resampling can be represented as a stochastic integral Jq{Z{X{s—), s) — 
X{s—))dYy(s). The evolution in species count due to both reaction dynamics and splitting is 

X^(i)=X^(0)+ V CY^^'i f\fiXA{s))ds)+ f F{XA{s))ds 

+ f iZiXA{s~),s)-XA{s-))dY^is); 
Jo 
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hence for the rescaled system = N ^Xa we have 



Xjv(t) - Xn{0) + C^^c^'f^^^'^f / iXN{s))a,N{l - XN{s)\Nd^s] + f FNiX{s))ds 

/^„^^=T V Jo J Jo 



(C,a,b)el 



+ / {N-^Z[NXm{s-),s) - XN{s~))d%{s) 



{N-^Z{NXm{s~),s) - Xn{s-))-i{NXm{s), N)ds 



Xn{Q) + MN^t) + / FN{X{s))ds 



+ / N-^{Z[NXn{s-),s) - NXN{s-))-i{NXN{s),N)ds 



(6) 



(6') 

Jo 

We still have 

Fn{x)^ N-+'-\^it{N){x)a^N{l-x\N, (7) 

(C,a,b)eli»" 

but now Mf^^j denotes the martingale formed by the second and fourth summand in (6) whose 
quadratic variation is 

[MN,^]t= Y N-^C^Y^-''(N''+'KfiN)f\xNis))a,N{l~XNis))b.Nds 



r-t 

+ / N-'^[Z(NXn{s-),s) - XXN{s-)f dY^{s). 
Jo 



(8) 



Note that since the two mechanisms are driven by independent Poisson processes, there is no 
quadratic covariation contribution. Since Ep^ . [Z{x, s) -x] =0 for all s > and x e {0, . . . , N}, 
the infinitesimal mean still satisfies 



ds 



E[Xw(s)|J-t] 



E[FN{XN{t))\^t] = N-^'^\KfiN){XN{t))a,Nil-XMit)) 



b,N, 



(9) 



on the other hand, the infinitesimal variance now satisfies 



ds 



B[{XNis)-E[XNis)]f\J't' 



= ^E[[Miv].lJ-t; 



= Y N^^'"\^i^t{N){XNis))a,N{l~XN{s))b,N+'y{NXr,is),N)N~''a%{NXN{s)) (10) 

(C,a,())eXbaluxbia 



2.3 Possible Qualitative Behaviours 



In order to determine the role that the rate of the splitting/resampling mechanism may play, we 
first establish the possible behaviour of the system when N is large. The decisive quantity for the 
qualitative behaviour of the system is 



where 



SA lim SAiN), e^(iV) := 



c„2{N):= Y N''+''-'^K,f{N)+ sup "fiNx,N)N-^alf{Nx) 



xe[o,i] 



(11) 



(12) 
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and 



c^(iV):= N-^'^^fiN); (13) 

ea relates the magnitude of the variance due to the splitting mechanism (or possibly a faster set 
of balanced reactions) to the magnitude of the drift due to reaction dynamics. If they are of the 
same order of magnitude, then the rescaled process will converge to a diffusion. In other words, if 
Ea € (0,00) then we can assume (by rescaling time as necessary) that both scaling constants (12) 
and (13) satisfy Co-2 = lim^v-j-oo Co-2 ( A^) € {0,oo),c^ = limN^oo Cfj.iN) G (0,oo), and Cg.2 = e^ic^. If 
Assumption 71 is satisfied, the noise of the splitting mechanism is such that the limiting behaviour 
of the system is diffusive, instead of being deterministic, as in (5) when only reactions are present. 

Proposition 2.1. IfsA e (0, 00), AssumptionT. holds for Xa = NXn, andXiq{Q) X{0) e [0,1], 
then Xpj ^ X as N ^ 00 in distribution on the Skorokhod space of cadlag paths on [0, 1], where X 
is a diffusion with drift and diffusion coefficients given by 

4>{x)= J2 C~^c>^{i~^}', H^)^ E C'~^t2x^{i-x)' + ^'a\x) (14) 

(a,b,C)e2:'"<' (a.bX)eX'"''' 

where for each (a, b, £ X*"" 

kf„ = lim N''+'"^Kf(N), 

for each (a, b, C) G l''"' 

kf^, = hm N''+''-\f{N), 

and for some 7^ > 

^^a^{x)= lim -f{Nx,N)N-^a%{Nx). 
If all reaction rates have standard scaling k!^ = k'^' N^~^°'~^''\ then k!^^2 = and a{x) = 7^(7^ (x). 

Proof. This is a direct consequence of standard theorems for convergence of Markov processes to a 
diffusion (see e.g. [ ] Sec 8.7.) based on locally uniform convergence of the infinitesimal mean and 
variance to the limiting drift and diffusion coefficients respectively, and convergence of jumps so that 
they disappear in the limit. Recall that the infinitesimal mean of the rescaled process Xn from (6') 
is given by (9) and its infinitesimal variance by (10). Since the process takes values in [0, 1], we can 
check convergence uniformly on the whole space, and moreover MN^j{t) = XN{t) — E[Xjv(i)], whose 
quadratic variation is given in (8), is then a square integrable martingale. For the contributions 
by the splitting mechanism, the convergence of the infinitesimal mean and variance, as well as the 
control of the jumps, are easy to check from the three requirements on the splitting mechanism made 
in Assumptions 6. and T. For the contributions by the reaction dynamics the convergence of the 
infinitesimal mean and variance, and the control over jumps, follow from the scaling properties of the 
counting processes used in their representation and from the fact that the rates for these counting 
processes are Lipschitz and bounded. These same conditions have been checked, in the case when 
reaction rates have a more general form, for law of large numbers and central limit theorem results 
for rescaled population-dependent Markov processes ['H]. Alternatively, one could also check that 
the Markov process X^ satisfies all the conditions required for convergence of more general Markov 
jump processes to a diffusion as stated in Theorem 2.11 of [ ] and Theorem 3.1 of [ ;(>]. The only 
thing left to check is whether a diffusion with coefficients as given exists and is unique in law. 
This follows easily from the fact that the contributions to a{x) and 4>i^) from reaction rates are 
polynomial, and we have assumed that ct^(x) is Lipschitz. □ 

A diffusion may or may not hit its boundary points, but it never spends a disproportionate amount 
of time at any point in its range, including the boundaries, unless they are absorbing. Hence, we 
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really need to consider the behaviour of the process when either — > or e^i — > oo (as a function of 
an additional asymptotic parameter which will be discussed below in Section 3.1). The only remark 
we make when ea remains bounded away from and oo is that the behaviour of X at the boundary 
{0, 1} depends on the form for the limiting variance of the splitting mechanism. As a consequence 
of Assumption 2., and of the properties of the splitting variance at {0, 1}, we are only guaranteed 
that 0(0) > 0, < and a(0) = a(l) = 0. Hence, {0, 1} are neither absorbing nor natural, but 
it remains to determine whether they are entrance or regular boundary points. Further conditions 
on the reaction and splitting mechanisms for reaching the boundary (i.e. for {0, 1} to be regular 
boundary points) are guaranteed by interpreting Feller's test for explosion (see e.g. [8], Sec 6.2., or 
[14], Sec 15.6.). 

The diffusive case sa G (0,oo) separates two other types of behaviour. When « and ea ~ oo, 
the rate of splitting is either slower or faster, respectively, than prescribed by Assumption 7? Both 
cases lead to behaviour which exhibits a type of stochastic bistability, in which the system spends 
almost all of its time at two points, or very near them. This bistability is, in the two cases Sa ~ 
and Sa ~ oo, caused by completely different effects of the two stochastic mechanisms in our model, 
which we investigate separately in the next two Sections. 



3 Bistable Behaviour from Slow Splitting 

Let us consider the case sa ~ 0, and assume that time has been rescaled so that = limTv-foo c^{N) <E 
(0, oo) and Co-2 = limjv_>.oo c^2[N) « 0. In modeling this is a relatively conventional scaling, in which 
a small amount of noise (from balanced reactions and splitting) will affect the predominantly deter- 
ministic behaviour due to drift (of biased reactions) . A precise statement of this depends on how 
fast EAiN) = '^^^1^^^ approaches as a function of N, and we examine it more carefully by first 
introducing a separate perturbation parameter e and then relating it to the scaling parameter N. 



3.1 Small Diffusive Noise Effects 



The simplest way to model small diffusive effects is with an enforced separation of time-scales between 
reactions and splitting using a perturbation parameter. Suppose all the reaction constants k'^^{N) 
depend only on the scale of the system N and have the standard scaling k"'' = f^ab]\j-i-{a+b) 
some constants k^^. Suppose the splitting rate, in addition to N, also depends on a small parameter 
e > 0, so that the splitting rate is '-f{x,N,e) — e'^j{x,N) where 'y{x,N) satisfies Assumption 7?. 
The fact that the splitting rate is slower than diffusive is expressed in terms of the fact that we will 
consider the behaviour of the system as £ — >■ 0. In this case the quantitiy Sa defined in (11) is just 
a constant multiple of 

Ea '■= hm — 



where c^2 = 7^ sup^gjQ ct^(x) and = Y.{a. 



z.ab 



We could also assume the rates of balanced reactions depend on the additional parameter , in the 
sense that nf = e^Kf'N'^-^''+^^ for (a, 6, C) S T'^'^y In this case 

5^2 = ~^ti^+f sup a'^{x). 

However, if we make no special separation in the way balanced and biased reactions are scaled. 
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then the assumption of standard scahng k"'' — k'^''N^ {a+b) jnipUes that this is only possible if the 
parameter s satisfies = N~^, on which we remark further in the next subsection. 

By Proposition 2.1, for any fixed e > 0, the process obtained in the limit Xj^ is a diffusion 

with coefficients <t>{x) as in (14) and d^{x) = e^7^(T^(a;) (we will use the subscript e in the notation 
of the limiting diffusion to stress its dependence on the small parameter e). X^. is a solution of the 
stochastic differential equation 

dX,{t) = ^{X,{t))dt + e'ya{X,{t))dB{t), X, e [0, 1] (15) 

where i? is a standard Brownian motion, a classical case of a diffusion with small diffusion coefficient. 

For many such diffusions £ ~ will have little qualitative effect relative to e = 0; however, suppose 
that (f) has two stable and one unstable equilibria, and thus the potential $ defined by $ = — / is 
a double-well potential. Since is a polynomial, this is an assumption on the number and type of 
zeros of 0. Explicitly, we will assume that: 



3 < xi < X2 < < 1 : (}){xi) = 0, i = 1, 2, 3 and 4)' [xi) < 0, (t>'{x2) > 0, (p'ixj) < 0. 



(16) 



Recall also that Assumption 2. implies that at the boundaries we have 0(0) > 0, 0(1) < 0. As 
a consequence, is a process whose mean behaviour involves monotone convergence to one of 
two stable equilibria (determined by the initial conditions), but where the small amount of noise 
allows the process to switch from one equilibrium to the other, creating a bistable system. Precise 
statements of this behaviour are described by Freidlin-Wentzell theory for random perturbations of 
dynamical systems by diffusive noise, [12], which can also cover processes with metastability, [ ]. 
We will follow closely the notation of [13], as their results apply most directly to X^. We first need 
a transformation to handle the state dependence (T^(x) of the diffusion coefficient, easily done using 
see [7] Sec 5.6., or [21] Sec 2.5. 

For Xg satisfying (15), large deviation theory for Gaussian perturbations of dynamical systems, [""] 
Thm 5.6.7 and Exercise 5.6.25, says that deviations of X^ away from an e-sized neighbourhood of 
xi and X3 are characterized by the large deviation rate function for Xg given by the quasipotential 
(with respect to Xi and X2) 



^i,X2(0,7^) 



:— inf inf 

s>Q 5 



eeCi([0,s]), e(0)-x„e(s)=X2 , z = l,3 



where L is the action functional 



Hi, a 



This identifies the most likely paths which leave a neighbourhood of Xi or x^, since every path 
between xi and X3 of the one-dimensional Xg has to pass through X2- We can write L(^,^') in this 
form for all such paths because Xg is non-singular away from the boundaries, that is, a'^{x) > c, Vx € 
[j;i, 0:3] for some c > 0. If the diffusion coefficient were constant ja = 1, then L(^, ^') = (^' — 4>{^))^ 
and Ixi,x2 (0j 1) would be simply a constant multiple of the potential, Ixi,x2 {4>, 1) = '^{'^{x2) — $(0;^)), 
for i = 1,3. The quasipotential would be determined by the height of the potential barrier which 
Xf, needs to overcome in order to pass from one equilibrium to the basin of attraction of the other. 

To solve the variational problem in our case, we can use a transformation of the path space ^ = g{ip) 
to get an action functional of the form i(^, ^') = (i/)' — 0(V'))^, from which we can deduce the explicit 
form of the rate function Ixi,x2 foi' state-dependent 'ya{x). For any monotone function g which 
for all s is surjective from C^{[Q, s\) to C"'^([0, s]), we have 



lxi,x2{4>,icr) = inf inf 

s>0 1/1 



L{g{i^{u)),[g{i'{n))]')du 



V e Ci([0,s]), V(0) = g-\xi),i^{s) = g-\x2) 
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We take g which satisfies the (autonomous) first-order ODE g'{y) — icr{g{y)), so that 

Note that ^d{x) > 0,Va; e (0, 1) ensures that g is in fact strictly increasing on (0, 1). Let h{x) = 
g~^{x). Then, if cj) is the vector field of a double-well potential, so is a defined as a = ^^^j 
for the following reasons. Let yi = g^^{xi) ~ h{xi); these will be the equilibria for a, since 
a{y,) = (i>{g(yi))/l^{g{yi)) = 4>{xi)/^a{xi) = 0. As for their stability, we have 

,/ ^ 4>'{9{Vi))9'{yi)lCF{yi)- 4>{g{y.i))^d'{g{yi))g'{yi) 4>' {g{yi))g' (yi) 7,/ . ^^ 

{g\y^)) 7CT(g(yj)) 

where the first equality holds since (t>{g{yi)) — 0, and the second by definition of g. Therefore, for 
each i, the stability of Xi under the vector field is the same as that of yi with vector field a; we 
may therefore define A = — J a to be the (double-well) potential associated with a. Since L(^,^') 
is now in the form L{g{tlj), [giip]]') — — a{g{ip)))'^ , we can conclude that 

h,,.A4>,l^)=Iy,,yAaA) = 2{A{y2)-Aiy^,)) i = l,3. (17) 



We can now interpret the results of [ ] to characterize the behaviour of the process (defined in 
(15)) as £ — > 0. Let Di denote basins of attraction for the deterministic process (5) driven by the drift 
(j), that is, Di — [0, X2) 9 xi, £^2 = {2^2}, D3 = {x2, 1] 9 2:3, and Bc{xi) denote closed balls of radius 
c > around Xi,Xz such that Bc{xi) C Di,Bc.{x3) C -D3. If the wells of the transformed potential 
A are not at equal depth A{yi) 7^ ^(j/s), we will without loss of generality assume A{yi) < ^(2/3). 
Let 

T, = inf{i > : Xeit) e B,{xi)}, f, - inf{t > T, : X^) e B^ix^)} 

denote the first hitting time of the neighbourhood of the stable equilibrium with the deeper basin, 
and the subsequent first hitting time of the neighbourhood of the other stable equilibrium. Let 
be the time-scale on which transitions from 1)3 to the neighbourhood of xi happen, defined by 
P[Te > f3^\Xir{Q) = X3] = e^^, and the one on which the reverse transition happen, defined 
by P[Ti, > /3^\X^{0) = xi] = e~^. The next result establishes that the transition from one stable 
equilibrium to the other happens on a time-scale of order 0(e^ ^(My2)-A{yi)) -j with i = 3 and i = 1 
respectively, and that in the limit as e the transition times have an exponential distribution. 

Proposition 3.1. If (j) satisfies (10), then the transitions of X^ from D3 to Bc{xi) and from Di to 
Ba{xz) satisfy 

[i) limP[Te > tl3^\Xe{Q) ^xeD3\ = er^yt > 0, limP[fe > t/3^\X,{0) ^ x e Di] = e"*,V< > 
{ii) lim e^ln/Se = ^3,,;, (0, 7a) = 2(^(2/2) - ^(2/3)), lim e^ln/S^ = 4,,,;, (0, 7a) = 2(^(2/2) - ^(yi))- 



Proof, (i) is a restatement of Theorem 1 in [ ]. (ii) follows from Theorem 4.2 of Ch 4 in [ /], 
which states, that for any S > 0, hm£_j.o P[|£^ hiT^ — /^g^j;^ (^, 7(t)| > 5\X,;{0) = X3] = and 
limE_).o P[|e^ InTe — Ixi,x2{4'ilS')\ > 6\Xe{0) — xi] = 0, and from our explicit calculation of the 
value of the quasipotential in (17). □ 



The following result characterizes the long-term behaviour on the natural time-scale (determined 
by /3e) for transition to the stable point with the deeper basin. Let = ° for some a G 
(0, 2(^(2/2) - Aijj^))), so that i?E -s> 00 while Re/Pe -> as e -> 0. Again following [13], define the 
measure-valued process {vl)t>Q by 

^KI) = 1^ / f{Xe(.s))ds 
Re J8,t 
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for any (bounded) continuous function / on [0,1]. The measure approximates the law for the 
location of X^{T) on the time-scale T — f5^t. 

Note that if A{yi) < ^(^3) then the results of (ii) imply that Pe/Pe — > 00 as e — 0, so that 
mixeDi P[21r//3e > t\ ^^(0) = x] — )■ 1. Hence, in this case metastability is characterized by the fact 
that the transitions into the deeper well are on an exponentially faster time-scale, relative to which 
the transitions back into the less deep well will not be noticed. Let Pa;[-] denote P[ • |Xe(0) = x]. 

Proposition 3.2. For each x G I?3, continuous function f on [0, 1], and S > we have 



lim 

6^0 



sup Wt{f)-f{x,)\>6 



= 0, hm P 



sup \iytif)-f{xi)\>6 



= 0. 



Moreover, we have convergence in law on the space of cadlag paths (with the Skorokhod topology) of 
(t'|)t>o to a jump process {vt)t>o such that 

(i) (Metastability) If A[xi) < A{x^), then {vt)t>{) is given by 



Vt 



<5,,, t<T 
t>T 



where T is an exponential mean 1 random variable, 
(ii) (Bistability) If A{xi) — ^(2:3) and a sequence of transition times is defined by: — 0, and 
Tl = inf{t > T-^ : X,{t) e B,{xi)}, fl = inf{t > Tl : X,{t) E 5,(2:3)}, i = 1, 2, . . . 
then (ft)t>o is given by 

Sx3, T2i <t< T2i+l 



Vt = { ^' - - - -^i+i 1 = 0,1,2,... 

\ S'xi , T2i+l < t < T2i+2 ' ' ' 

where Tq — 0, and {Ti}iyQ are arrival times in a rate 1 Poisson process. 

Proof, (i) is simply a restatement of the main result Theorem 2 in [1.3], and (ii) is an easy extension 
of this result. Since A{yi) = A^y^), we have = /^e and the transitions from one stable equilibrium 
to the other happen on the same exponential time-scale. By Proposition 3.1 («), on the time-scale 
T — f3^t, in the limit as £ — >■ 0, T^ is exponentially distributed with parameter 1, and X'^{T^) E Di. 
By the strong Markov property of X^, the time increment to the subsequent transition T^ — T^ is 
independent of , and the same Theorem implies that on the time-scale T = P^t = P^t, in the limit 
as £ — >■ 0, T^ — T^ is also exponentially distributed with parameter 1, and X'^{T[) E D^. The rest 
now follows from the same arguments as in the proof of Theorem 2 in [ .;]. □ 



3.2 Finite-System-Size Effects 



The above results relied on using an additional parameter s to separate the scaling of the noise 
from the scaling of the drift, obtaining a diffusion approximation for the limiting process first, then 
applying large deviation techniques for the diffusion (15) with small perturbation coefficient e. A 
priori, there is no reason why the limits need be taken in that order. Another approach is to apply 
large deviations techniques directly to the rescaled process Xn = N~^Xa, and obtain results that 
describe the large time-scale behaviour of X^r relative to the equilibrium points of the limiting drift 
(16). It is natural to compare these results to those for the associated diffusion with small diffusion 
coefficient. We will identify the exact relation of time-scales of the reaction system and the splitting 
mechanism for which large deviation rates of these two methods can be compared. 
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This question is most easily answered when the reactions and the sphtting/resamphng mechanism 
make only unit net changes at each step, so that Xa is a birth-death process. Assume, as before, all 
the reaction constants have the standard scaling k'q' — K°^iV^~(°+''), and assume again the splitting 
rate is of the form 7(x,iV, e) = e^7(x,iV), except now the parameter depends on N as well. 
Since, by Assumption 6., the change due to splitting is unbiased, we have Px,x+i — Px.x-i — 5, and 
the splitting variance is a% = 1. As earlier 7(x, N) is assumed to satisfy condition in Assumption 7!, 
that is, sup^ \jix,N)N-^ - 'y^a^i§)\ 0. 



Suppose Xa is a Markov jump process with rates Nf^{x)dt = P[Ar^(i + dt) = x + \\XA{t) — x\ 
and Nf^{x)dt = P[X^(t + dt) = x — l\XA{t) ~ x] such that lnf+,lnf_ are bounded Lipschitz 
continuous functions, and Xpf = N^^Xa is its rescaled version. Then, according to the Freidlin- 
Wentzell large deviation theory for Markov jump processes, [JV] Theorem 6.17, since transitions 
between the two stable equilibria xi,X3 of are uniquely achieved by crossing the potential barrier 
at X2, the deviations of Xn away from neighbourhoods of xi and X3 are characterized by the large 
deviation rate function for Xjv given by the quasipotential (with respect to Xi and X2), 



inf inf • 

T>0 5 



mu),e{n))du 



^eC\[0,T]), ^{0)^x,,aT)=xA, z = l,3 



where £ is the action functional in variational form 

£ix, y) = sup {ey - {f+{x){e' - 1) + r_ [x)[e-' - 1)) } 


determined from the jump rates of the process f_|_ and f_. Calculus of variations results, see 
Theorem 11.15, give an explicit expression for the quasipotential as 



In 



f_(a;) 
r+{x) 



dx. 



1,3. 



(18) 



If Xa is a birth-death process whose rates r+{x),r^{x) are such that r^{x) = N~^r^{Nx) f^+(a;) 
and r^(x) = N~^r-{Nx) — >■ f-{x) uniformly in x G [0, 1], then the logarithmic moment-generating 
function gN{x, 0) of the jump measure ijln{x, ■) — r^{x)6\ + r!^{x)S-i, for fixed 6, also converges 
uniformly in a; G [0, 1] 

57v(x, 0) = J (e«- - l)fiN{x, dz) = r'l {x){e' - 1) + {x){e-' ~ 1) 

r+{x){e'> -\)+r^{x){e-'' -\)= \ {e^^ - \)^L{x,dz) = g{x,Q) 

N^OD J 

to the logarithmic moment generating function of the jump measure /i(x, •) = r^{x)di + f^{x)6^i. 
Since the Legendre transform i]\[{x,y) of gN{x,y) has the explicit form 

iN{x,y) =sup{9y - gN{x,9)} 
s 

-y + ^ 4r^(a::)r^(a;) 
2r^{x) 



for fixed y, we also have uniform convergence in a; G [(5, 1 — (5], for any S > 0, 

^N{x,y) = sup {0y - gN{x, 6)} — > sup {Oy - g{x,e)} ^ e{x,y). 

9 N^oo 

Consequently, the large deviation behaviour for X^ — N^^Xa is determined by the same action 
functional £(x, y) and exit times in terms of the same quasipotential ixi,x2 {f+i ^^_), z = 1, 3 as above. 
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For the system of reactions and splitting, birth and death rates for the process Xjv, r+ and r_ 
respectively, are of the form 



(a,6,l)el 



a.N \ ^' / b.N 



We wish to obtain results for the time-scale of exit from a neighbourhood of a stable equilibirum 
for the rescaled process Xn that are analogous to those for obtained in Proposition 3.1. To this 
end, we will have to make some assumptions about the behaviour of and r_ in order to use the 
quasipotential la;;^^;^ (f+, r_). Let ^'^d Pen denote time-scales of the transitions of the process 
Xn from D3 to Bc{xi), and from Di to Bcix^) respectively, in the analogous way as /?£ and were 
for the singularly perturbed diffusion. The next result establishes the time-scale of transition for 
Xn from one stable equilibrium to the other. 

Proposition 3.3. // Xa is a birth-death chain, whose rates satisfy: 

^^^ffl^f-(-) uniformly tn [0,1] (21) 

such that (f> ~ r^ — satisfies (16), then the mean times jS^^ and for transitions of Xn from 
to Bc{xi), and from Di to Bc{x^) respectively, are given in terms of ij^.^j;^ (f_|_, f_) from (18) by 

lim lim 



Proof. This is just the statement of results for the exit problem for the jump Markov chain Xn in 
terms of its quasipotential, obtained by Freidlin and Wentzell [12]: see Theorem 1.2 and Theorem 
2.1 of Ch 5, the discussion at the beginning of Section 4 and Theorem 4.3 of Ch 5; also see Theorem 
5.7.11 of Ch 5 in [ ]. Uniform convergence of the action potential, that is the Legendre transform 
£n, is necessary in order to express the quasipotential ixi,x2 in terms of the limiting rates f+,f_. 
All of the assumptions on the equilibrium points of 4>{x) = r+(x) — f_(x) in (16) are also necessary, 
since (j) determines the fluid limit of the jump Markov chain Xj^. □ 



Finally, we can establish the time-scale separation under which the switching results for the rescaled 
jump process X^ and the diffusion with the small diffusion coefficient can be compared. 

Theorem 3.1. If the reaction system has increments of size {1,-1} only, its rates have standard 
scaling k'q' — k'^' N^~'^°-^^\ its limiting drift <f> satisfies (16); and, if the splitting mechanism has 
increments of size {1,-1}, its rate is £^7(2;, iV) where j{x,N) satisfies the Assumption 7* and 

Ne% 1; 

then, results based on large deviations for Xj^ in Proposition 3.3 are more informative than re- 
sults based on large deviations for the diffusion X^ with the small perturbation parameter ejy in 
Proposition 3.1, that is 

Proof For r^{x) = N~^r+{Nx), and r^{x) = N-^r^{Nx) by (19)-(20) we have 
r+i^)^ J2 ~<{x)a.Nil-x)b,N + \N-'e%j{Nx,N), 

{a,b.,l)el 

r^{x)= J2 ii-\ix)a,N{l-x)b^N + lN-'s%jiNx,N). 

{a,b.-l}el 
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Since ^{x, N) is such that \^{Nx, N)N ^ — 7^(7^(a;)| — > uniformly in x e {0, . . . , 1}, then given 
that Nef^ -> 1, we have uniform convergence of — >■ f+ and — > f_|_ to 

(a, 6, 1)61 

= ^ k'^_\x^{l ~ xf + J7'^'(a;) 



Let c^(a;) = l-H^ 



SO 



f+(a;) - 4){x) 
ui{x) = — 



a;(a;) f+(a;) — f_(a;) ^(2;) 



f^{x) E(aA-i)ex^-''i^°(l-^)'' + 57''^'(2^) 

and (18) imphes that ixi.x2 {^+: = J^^ In (l — aj(x))(ia; satisfies 

(f){x)dx _ /"^^ (j){x)dx 



J2 k'i\x^il~x)'' + ^J^^^x)-'"'•"'^''^''^^- ^ E Aifa;'^(l -a;)'' + 172^2(2;) 
(a,f),-i)ex (a,fca)ei 

On the other hand by (17) and the fact that g'{y) — jS-{g{y)) we also have 

4...(0,^)=-2r«(y)rf.=-2r^(44^ > ^_.(f„f_). 

4. -'a. 7c^(.9(y)) Ja;. ^T<^^{x) 

Hence, if -/Ve^ — ^ 1, we get a comparison using quasipotentials for Xjv and of the time-scales for 
transitions between stable equilibria, as 

□ 



If £^ = N~^, transitions between stable equilibria are more often due to finite-system-size effects 
than due to the effects of an additional mechanism. This is understandable in light of the fact that 
the diffusion is a limit of the rescaled process in which the contribution of any subdiffusive 
noise disappears. As remarked earlier, when — iV^^, we could use this informally prior to 
obtaining a diffusion limit . If, for rates of balanced reactions we write k^''(7V) = N^-('^+'>)f;^f> = 
£2jY2-(a-i-6) -ab^ ^T^^^ ^j^g (jiffi^gion cocfRcicnt would become a^{x) — Z](a 6 c)e/'"'' '^C m^"'-"'^ ^ 
x^ + 72 0-2 (a;)). However, even this 'adjusted' diffusion coefficient would not change the conclusion 
of Theorem 3.1, since the contribution of the rates from biased reactions is still missing from the 
quasipotential of Xg. 

If <C £^ <C 1, it is clear from Theorem 3.1 that the noise of the splitting is the dominant factor 
in effecting the transitions, while if ^ N~^, the noise from reactions dominates, and both rates 
f+,f_ and the quasipotential ixi,x2 determined by the reaction system only. 
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3.3 Example: Bistable Behaviour from Slow Splitting 



Here is an example of a simple reaction system that yields a limiting system with a double-well 
potential. 

A B (22) 

b"X A (23) 

A + B -4 2B (24) 

iA + B^A 'iA (25) 

The trimolecular reaction (25) produces a term in the drift which is cubic in XA^ which is needed in 
order to obtain the three desired equilibria. With standard mass-action scaling = A^^~(''"'"^)k^^, 
the limit of FN{XN{t)) = E[XAr(t)] = 'E[XAit)/N] G [0, 1] as iV ^ cx) is 

4){x)= lim i^Ar(a;) = + - 2^) - K"ia;(l -x) + Kf a;2(l -x), a;e[0,l]. 

N^OO 

With the special choice of k^-^ — k^^ — 1, k^^ — kf^ = ^ we have 

^{x) = ^ (3 - 22x + 48x2 - 32x=') = -f " ^) " ^ " I) (26) 

with two stable points at xi = ^ and X3 = | and one unstable point sX X2 = \ for the system, 
and thus $ = — J (/) is a double-well potential. Since (j) is antisymmetric about the line x = ^ the 
potential can be expressed as 

$(x) = i(2x-l)4-^(2x-l)2 + C, 

which is symmetric about the line x = ^ and thus $ has equally deep wells '^{\) — *i>(|). 

This system bears resemblance to the so-called Schlogl model [ ] , which consists of four reactions 
A + 2X ^ 3X, B ^ X, with the resulting drift for X cubic. In [ ] the authors formulate the 
Kolmogorov forward equation (chemical master equation) to analyze the stochastic model for this 
reaction system. 

For this example we take the simplest split ting/ resampling mechanism (Example (Bern) in Sec 2.2) 
in which at each split an error in the molecular count of A from the parent to the daughter cell is 
at most 1. Its rate is 7(x, N) ~ j{N)x/N{l — x/N) and its probabilities are Px,x+i = Px,x-i = 1/2 
for X 0,N, and po,o = Pn,n — 1- Note that here the factor ^{N) will depend on TV, but is state 
independent. This mechanism can also be represented in terms of reactions as 

^ + B^"^^^2A A + B 25 (27) 

We stress that this representation (27) of the resampling in terms of reactions is done merely to 
illustrate the mechanism in a similar way to the reactions, and is not to be confused with an actual 
set of biological reactions as in (22)-(25). This can be done in the particular case of Moran-type 
resampling, since the rates of this mechanism depend on the product of both the count of A and of 
B. This is a consequence of the fact that each resampling event picks either one molecule of A or 
one molecule of B with probabilities relative to their proportions in the cell, and replaces it in the 
daughter cell with a random choice of either A or B with equal probability. 

As shown in Sec 2.2, if we choose the splitting parameter to be j{N) = ^e^N"^ for some small 
constant > q, then 7(x, N) satisfies all the conditions of Assumption 7*, and the limiting process 
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Xe satisfies the stochastic difi'erential equation with drift (2G) and diffusion coefficient e^^x{l — x) 

dX,{t) = ^ (3 - 22X,{t) + iSX^it) - 32X^{tf^ dt + e\J X^{t){l~ X^{t))dB{t). (28) 

To find the value of the quasipotential for this problem we find the transformation of the potential 
via a{y) = 4>{g{y)) /a{g{y)), where g is the solution to g'{y) = a{g{y)) = ^/g{y){l ~ g{y)), given 
explicitly by 

.(y)=cos^Q(y-|))=cos^ (Fi 

We chose the constant of integration so that 5(0) = ^, and g{—y) = 1 — g{y)- The inverse of g is 
given by 





TT TT 




2' 2 



TT TT 

2^6' 



h{x) ^ g ^(x) = 2arctan 1^ - y - - ij + -, x € [0, 1] 
hence, the transformed equilibrium points yi — h{xi) are 

2/1 = 2arctan(-y3) + | = 2/3 = 2arctan ^ - 

and 

TT 

2/2 = 2arctan(— 1) + = 0- 

Note as well that the wells of the transformed potential are of equal depth, which follows from the 
fact that a is an odd function 

^ a{g{-y)) a{l - g{y)) d{g[y)) 

and thus A = — J a{y)dy is an even function. Since yi — — j/3, j/2 — 0, and A(yi) = Aly^) ~ 0, we 
have A{y2) = J^^^ a{y)dy with a rather complicated expression 

0(cos^(|-|)) 1 /-^/^ 3- 22^ + 48x^-32x3 

A{y2) = / —, ^-77^ ^dy = - / r dx = 0.0913 

Jo cr(cos2 (f - f)) 3 7i/2 x(l-x) 



By Proposition 3.1 on a time-scale of 0(e"^ 2A(y2))^ ^^e process exists a neighbourhood of the stable 
equilibria xi — \,X3 = |. Symmetry of A around j/2 = implies that we are in the bistable case 
(ii) of Proposition 3.2, and the occupation measure of the process X^ converges to the occupation 
measure of a two-state Markov chain, which transitions between states {3;, |} with equal rates. 
Figure 1 shows an exact simulation of a sample path of the rescaled process Xjv = N~^Xa with 
choice of parameters N — 1500, — 0.02; since ^ 1/N, we expect the e-perturbation of the 
limiting diffusion to be driving the switching. Indeed, the process appears to be spending most of 
its time in neighbourhoods -60.1(2^1) U Sq. 1(^3)1 switching between them at the approximate rate 
R = e-^"'2A(a2) ^ 0.0001083. 

If we take <C l/N, then transitions between stable equilibria are based only on the scaled rates 
for the reaction system (22)-(25), 

r^[x) — 1 ~ X + —x''(l — xj and r-[x) — x + —x(l — xj. 
By Proposition 3.3 the values of the quasipotential for the birth death Markov process are 

/f (x)\ n ( x^^-xiX-x) \ 

In ]dx^ In — ' ^ 0.006713, 

\r+{x)J Ji \l-x+fx^il-x)J 
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500 1000 1500 2000 0.25 0.75 



Figure 1: Sample path Xi\[{t) (left: x-axis=t, y-axis=XN{t) — N^^X{t)) and its occupation density 
(right: x-axis=state space of C [0,1], y-axis= proportion of time Xjv spends in each state 
by time t = 2500) from the system (22-25) with birth-death splitting, under standard mass-action 
scaling for reactions and 7(5, TV) = ^e^TV^ (parameters N — 1500, £^ = 0.02). Dashed red lines 
indicate quasi-equilibria at 1/4 and 3/4. 




Figure 2: Sample path XN{t) (left: x-axis=t, y-axis=Xjv(i)) and its occupation density (right: 
x-axis=state space of Xn C [0,1], y-axis= proportion of time Xn spends in each state by time 
t = 4000) for the system (22-25) with birth-death splitting, under standard mass-action scaling for 
reactions and j{e,N) = ^e'^N^ (parameters N — 500, = 2 x 10""*). Dashed red lines indicate 
quasi-equilibria at 1/4 and 3/4 as above. 
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and 

[ In ( ^^4^ I dx = 0.005534 

Note that here the values for the quasipotential are no longer equal, and the process will take longer 
to get out of the neighbourhood of the equilibrium xi = j. Figure 2 shows a simulation of a sample 
path of the rescaled process Xn for 7(7V) = e'^N with the choice of parameters N = 500, but 
= 2 X 10""*. In this case 1/A'' and we expect the transitions to be due to noise from the 

reactions arising from finite-iV effects. Based on the above calculation we expect the process to be 
switching away from Bq i{xi) at rate R = e~^*^i'^2 = g-3. 356629 _ 0.035 and away from -Bq.iI^^s) 
at a rate R' = e~''^*=^3'"=2 = g-2. 769957 _ 0.062; indeed, the time spent near xi is appreciably larger 
than the time spent near 2:3. 

We make a particular note that the reaction system considered here is very sensitive to the exact 
values given for the reaction constants; a small change in these would preserve the double-well 
potential, but would lead to nonequal depth of the two wells for the quasipotential, and hence 
instead of a limiting bistable behaviour would lead to a to limiting metastable behaviour as in case 
(i) of Proposition 3.2. In the next section we discuss the conditions on the scaling of the reaction 
and splitting/resampling which yield behaviour that can also be described as bistable, but where 
the underlying mechanism is qualitatively different and the restrictions on the reaction system are 
negligible. 



4 Bistable Behaviour from Fast Splitting 

We next consider the case ea ~ 00, and assume that time has been rescaled so that c„2 = 
lmij\j^ao Ca^{N) S (0,00) and — limjv-i-oo c^(-^) ~ 0. This is a more unconventional scaling, 
in which the noise (from balanced reactions and splitting) overwhelms the contribution due to the 
drift (from biased reactions). 

One way to model this with a diffusion would be to introduce a time-scale separation with an 
additional small parameter e in the scaling of all reactions rather than in the rate of splitting. 
Suppose all reaction constants scale as — eK'^^N^~'^°'^''\ while the rate of splitting 'y{x,N) 
satisfies Assumption T. . For any fixed e > 0, the resulting limit of the rescaled process would 
be 

dX%t) = e4>{X%t))dt + ^d{X^t))dB{t), X" £ [0, 1] (29) 

where _B is a standard Brownian motion, and we have the case of a diffusion with a small drift. Note 
that although o-^(O) — (t^(1) = (by Assumption 7.), the boundaries {0, 1} are not absorbing, since 
there is at least one biased reaction that allows escape from either boundary 0(0) > 0, </)(!) < (by 
Assumption 2.). Other than at the boundaries the contribution of the drift is essentially negligible, 
and X^ is approximately a martingale. Most attempts to escape a boundary are followed by the 
return to the same boundary point, only some end up at the opposite one. In the limit as e — > 0, 
the rate of escapes from the boundaries for X^ vanishes, and there is no switching. 

However, under the right conditions, the limit of the original rescaled process will spend almost all 
of its time at one boundary or the other, switching between the two on a reasonable time-scale, 
creating again a bistable system. How the effect of the attempts to escape the boundary appears in 
the limit depends on the rate of the attempts, and the time spent between the boundaries. In order 
to make a precise statement we need to examine the behaviour of the rescaled process X]\[ directly 
and specify a general set of conditions for a Markov jump process to exhibit this type of switching 
behaviour. 
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4.1 Stochastic Switching 



The unsealed process Xa is a Markov chain on {0,1,..., N} with transitions that are due to the re- 
actions (a, fe, C) G 2^, as well as the splitting mechanism with distribution Px,y, {x, y) £ {0, 1, . . . , A^}^. 
The rates of these transitions from Xa = x are equal to ^ q^j^ A^''(a;) = ^ ^jg^; K^''(iV)(x)a(-/V— 
x)b from the reaction system, and 7(x, N) from the splitting, respectively. We denote the total com- 
bined rate of Xa from i e {0, . . . , n} to j € {0, . . . , n} by 

- ^ ,(7V)Za(iV - l), + 7(z, iV)p,,,. 

(a,fe,j — 

Transitions due to splitting can have jumps whose size can in principle be as large as — 1 (such 
as those of the Wright-Fisher splitting process example in Section 2.2), although with very small 
probability. However, a splitting mechanism is absorbing at {0, N}, po.o = Pn,n = 1, and the rates 
of jumps off the boundaries x S {0, iV} are created by reactions using only molecules of B (for 
a; = 0), or using only molecules of A (for x = N), with rates 

(o,fcj")ei (a,o,N-j)ei 

By Assumption 2. in Section 2.1, there exist j, j' G {1, . . . ,N — 1} such that roj, r^j' 7^ 0. The 
leading powers of N , max(o > and niaxj^ o,7V-j")ei{a} > respectively, will determine 

the rate at which attempts to counteract absorption at the boundaries happen, and in particular, 
this implies that rpj, r^vj — >■ 00 as A'^ — >■ 00 (allowing for upcoming condition (30)). 

Define an excursion of Xa to be any segment X^(i), t e [^1,^2) such that X^(<i— ), X^(t2) € {0, 
and XA{t) ^ {0,iV} for t G [^1,^2)- Call an excursion on [ti^t^) "successful" if X^(ti— ) 7^ XA{t2), 
and "unsuccessful" otherwise. For < j < A^, let Tj :— inf{t > : X^(i) = j} be the first hitting 
time of state j , and let tqjv = tq A tjv denote the first hitting time of either boundary state. Let 

ejo = E[ro,Ar|X^(0) = j,X^(to,7v) = 0], e^N = E[ro^Ar| Xa(0) = j, X^(to,7v) = A^] 

be the expected hitting time of the two boundaries from j, and ■Kjn be the probability that an 
excursion from j hits the A^ boundary first 

Ti,N = V[Xa{to,n) - A^l Xa{Q) = j], 

and thus tTjo = 1— tTjat is the probability it first hits the boundary. The values of {cj. , T^jN}j£{i....,N-i} 
can be determined by setting up and solving the appropriate linear functionals of the generator for 
the Markov process Xa\ explicit expressions, however, may be hard to come by for general processes. 

Excursions of Xa depend on transitions from both reactions and the splitting mechanism. However, 
if the noise overwhelms the drift, then at each step in the interior transition rates are dominated by 
those from the balanced reactions and the splitting mechanism. In particular, this will imply that 
in the interior Xa behaves approximately like a martingale, and will allow us to approximate the 
probability of switching from one boundary point to the other in terms of the relative rates of biased 
reactions versus balanced reactions and splitting. We will estimate Cjo, ej^ jT^jn in an example to 
come, and exhibit more explicit conditions than the ones below in the case when the reactions and 
splitting yield a birth-death process for Xa- 

We first state general conditions under which the rescaled process Xn = Xa/N can be approximated 
by a simple Markov jump process. Suppose that there exists two scaling parameters: the order of 
magnitude of the rate of reactions on the boundary wjv 00, and a time scaling parameter /3jv > 



22 



for the rescaled process Xjv, such that: 

^ H ^ ' 77" ^ ^- (3^^) 



/^AT 2_^roj TTj^ -> foi, /3jv 2^ ''JVj TTjO fio (31) 



XI '^OJ '^J^' m ''OJ ^JO' ^^^^ XI ^JO' X '^^J '^oN (32) 



with f+, f_, foi, fio G (0, cx)). Since roj^r^j — >■ oo, there is no need to change the time-scale for 
the process. These conditions imply that there are many excursions in any finite time interval [0,i], 
only a small fraction of which are successful, and during which the total time spent is very small. 
Consequently, the rescaled process will spend most of its time on one boundary until the first time 
a successful excursion takes it to the other boundary. Let = inf{< > : X^(t) = 0}, and 



Ti^^mf{t>f^-^:XAit)^N}, = inf{i > : X^(t) - 0}, i-1,2... 

be a sequence of times at which Xa first reaches a boundary different from the one where it was 
most recently. Also, define the measure-valued process {vP)t>a for some pM > such that ^ — )■ 

by 

^tU) = — / f{XN{s))ds 
Pn Jji^t 

for any (bounded continuous) function / on {0, . . . , 1}; this {v^) approximates the law of the 
location of the rescaled process X^it) = XA{t)/N on a short time interval of length p^. 

Proposition 4.1. If Xa satisfies (30)-(32), then 

lim P [TAr - f^-i > tl3N] = e-^«i', lim P [T^+' -f'^> tpi,] = e-''°\ Vt > 0, 

and we have convergence in law on the space of cadlag paths (with the Skorokhod topology) {ly^ jfX) 
(^'t)t>o to a jump process 

where {r2i+i —T2i}i>Q o,nd {121+2 —T2i+i}i>vi tf"e two independent sequences of i.i.d. exponential 
variables with rates foi and fio respectively. 



The rescaled process Xn can therefore be approximated by a jump Markov process (J(i))t>o on 
{0, 1} with transition rates fgi from — 1, and fio from 1 — > in the following sense: the occupation 
times of X^r on {0, 1} converge to the respective occupation times of J, and the times of successful 
excursions of Xn from 0^-1 and from 1 -> converge to the respective transitions of J. We cannot 
expect a stronger kind of convergence than stated, since for example, convergence in law of Xjv 
to J in the Skorokhod topology is precluded by the fact that for arbitrarily large N, there remain 
unsuccessful excursions of Xn that stray from their originating boundary by a distance which is 
bounded away from 0. 

A different set of conditions from those in (32) for the length of excursions away from the boundaries, 
where in the limit we get four non-zero limiting constants eoi, eoo, eio, en, would imply convergence 
to a limiting process which spends a nontrivial fraction of time away from the boundary. The limiting 
process would behave similarly to a diffusion with 'sticky' boundaries (see [ ] Sec 15.8 C). 



Proof of Prop. 4-1- For each i > 0, define a sequence of times after at which excursions from 
start a]^ and end f^* , by letting = T^, and for i' = 1, 2, . . . 

a]:;' = inf{f^''-i < t : XA{t) ^ 0, X^(t-) = 0}, f}^' = M{f;^''' < t : ^^(t) - 0, X^(t-) ^ 0} 
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and let s{i) — inf{i' > 1 : f^* > T^^} be the index of the first excursion from that is successful, 
hence f^"^*'' = T'}^^. Note that ^^(t) = 0,Vt G [f]^'' ) and that J2i'<sii)i'^N' ~ 

is the time spent at between successful excursions; while Xa 7^ for t G [a]^ ^f]^ ) and thus 
Si'<s(i)('^w' ~ ^^N ) ^^"^ time spent on unsuccessful excursions. 

Consider the time interval [T^,Tj^"'^] — Uj/<s(i) [ct^* ,f^* ) from which subintervals for unsuccessful 
excursions are excised. Excursions from are started at overall rate 7'oj', and since excursions 
whose first step is to j are successful with probability ■KjM^ successful excursions are started at rate 
Ej^OjTTjw So, 

Wn ■= ^pf'^'^ -Tl] - ^ (r^* -a'^)^ exponential f y^jQj-^jN 

i'<s{i) ^ 3 

and 

s(i) ~ geometric — 

Also, for any i' < s{i), the unsuccessful excursion times f]^* — a]^^ are independent and identically 
distributed with 



-^n]^2_, ^^^E[ro,Ar|X^(0) = j, X^(to,w) = 0], 



while T]^^ — (t]y*'''* is a subinterval for a successful excursion with 



Let 



C^:=E(#-^^') and 5;,:=T^i-a:v^«, 

i'<s{i) 



SO that — = VF^ + C/^ + 5)^. Assumption (31) implies W]^//3n => exponential(foi) as 

N ^ 00. We next show convergence for both Uj^ and S}^ in probability as — 00, which 
wiU imply that (T^^ - fi^)/^^ exponential(roi)- 

We first note that 

therefore, E[S']y//?Ar] — ^ 0, since the first fraction converges to and the second to 0, by (30) 

and (32) respectively. Similarly, for each unsuccessful excursion 1 < i' < s(i) 



and since s{i) is geometric. 



E[s(i)] = „ ' = UN 



We have 



E[C7^] = E 



(Tjv - ctn , 

i' <s{i) 



< E[s(j)]E[t^ - c^at j = ^ ^ • 
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and so E[[/^//3Ar] -> 0, since by (31) the denominator converges to fgi, and by (32) the numerator 
goes to 0. Hence for any (5 > we have P[Slj > 6] < -> and P[Wj^ > 6] < 0. 



A completely analogous proof shows that (T^ — T^)f3N exponential(fio), and the claim about 
the probability measure i^t is immediate from the fact that E[?7^ + S^^] — > 0. □ 

To verify condition (30) one only needs to use the rates of biased reactions on the boundary. For 
(31), note the fact that if not for biased reactions, the process would be a martingale; if the rates of 
the biased reactions are overpowered by those of the balanced reactions and splitting (as quantified 
in (31)), then the process is approximately a martingale. Conditions in (32) predominantly depend 
on how fast the rates of the balanced reactions and splitting are, as they determine the length of 
excursions of the process. 

These conditions are the easiest to verify when the reactions as well as splitting/resampling mech- 
anism make only unit net changes at each step, so that Xa is a birth-death process with r^- = if 
|i — j'l > 1. In this case one can specify more precise conditions on the rates rij that will ensure that 
(30)-(32) hold. We consider the case when the time is already rescaled, that is ^S^v = 1, and the rate 
of reactions on the boundaries is ujn = N . We use the following notation for birth and death rates: 

(^) 

with eAr(z) quantifying the strength of the bias at state i (we stress its dependence on N via transition 
rates r±{i)). 

Proposition 4.2. // Xa is a birth-death chain whose rates satisfy: 

^^f+e(0,(X3), !:^^f_e(o,c»); (33) 

N-l 

J2\eN{i)\^0; (34) 
1=1 

then the conditions (30)-(32) hold with ujn = N, f3]y — 1 and foi = ^lo = 



Analogous to the general case, (33) depends only on the rates of biased reactions on the boundaries, 
(34) reflects the fact that off of the boundaries the drift of the biased reactions is much weaker 
than the noise of the balanced reactions and splitting, and (35) is an estimate on the speed of the 
balanced reactions and splitting. 

Proof. (30) is immediate from (33) and ujn = N. To verify (31) we solve for TTjN,j S ■ ■ ■ , N —1}. 
Lemma 4.1. // (34) holds, then Ntti^ 1 and A^7r(^_i')o ^ 1- 

Proof. Let (p be such that lp{Xa) is a martingale, that is, let (p{x) = 'E[(p{Xa{to.n))\X{0) — x] 
fora;S{l,...,A^ — 1} and (p{0) — 0,(p{l) = 1. Standard result for birth-death processes, using a 
recursive equation for — vi^) ~ vi^ ^ l)i gives 



i=l i—1 j 
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By the Optional Stopping Theorem for the stopping time tq.at, 
so TT,N = Mi) - ^{0))/{^{N) - ^{0)) = ip{i)/^{N), and 



where c{N) = ^ Ef=i n;=i(l + ^nU)). 

The condition (34) imphes that sup]^<j<^_]^{|eAr(j)|} 0, so let A'o be such that VA^ > A'o and 

e{l,...,N~ 1}, |ejv(j)l < 1/3- Since Vx S [0, 1/3), 1 - a; > e-"=-^', and Vx G M, 1 + x < e"^, we 
have that uniformly for all 1 < a, 6 < iV — 1, where N > Nq 

b N-1 N-1 /n-1 \ 

i[(^+snU))< n(i+i^^wi)^ n^""^'^'-°^p Ei^^'O')! (36) 

j=a j = l j = l \i = ^ I 

and 

6 6 N-1 

11(1 + e^ii)) > 11(1 - |e^.(J)|) > n (1 - 

i=a J=a i=l 

Af-1 / N-1 \ 

> n exp (~|e^(j)| - kA.(j)n = exp - ^ (|e^(j)| + |e^(j)r) (37) 
hence A^ttuv = l/c(A^) ^ 1- 

To get -/V7r(jv-i)o ~^ 1? if we flip the state space by letting I ~ N — i, then the new boundaries are 
= iV and N — 0, and we get a birth-death process Xa whose rates are precisely the flip of those 
for Xa- That is, the rates of Xa are: f+(t) = r^{N — i), f_(i) = r+{N —- i), and their ratio is: 

giving the same product of ratios as for the original process. 

JV-l N-1 

^(l+^^^('^))- n(i+^^w)- 

1=1 3=1 

Hence, the exact argument above now applied to Xa gives -/Vttjjv ~ -^'"'(Ar-ijo — > 1 as well. □ 

Once we have the result of Lemma 4.1, it is immediate that Nttin — > 1, A'^7r(jv-i)o ™ply (31) with 
foi = f+,f 10 = f_. 

To verify (32) we next solve for ej,j G {1, . . . , — 1}, where Cj = E[ro,Ar|X^(0) — j] for j G 
{1, . . . , iV — 1}, and bq — cn — 0. 

Lemma 4.2. // (34) and (35) hold, then Nei and Ncn-i 0. 

Proof. The expected time of an excursion satisfies the recursion 

1 r_(i) r+{i) 



r- (i) + r_|_ (i) r_ (i) + r^. (i) r_ (z) + r+ (i) 
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which gives 

r+(i)(e,;+i - e^) - r_(i)(ei - e^-i) = -1; 
letting /(i) = — ej_i gives the recursive equation 

/(^ + 1) = — ^ + = — ^ + (1 + £A' 



Note that except for the — ^r^j^^ term this is reminiscent of the recursion for ^p{i) = — !)■ 

Hence, 

fc-i fe-i fc-i 

/(fc)-/(i)n(i+^^w)-Er7iy n (i+^A^c?))- 

To find /(I) = ei — eo = ei we impose the condition X^iLi /(*) = ejv — ep = and get 

^ TV fc-l k-l V I / ^ 

n (i+^A^o-)) / En(i+--o-)) 

Let ?7Ar = supj<^ h<jv-i I X^J^ + ^nU)) ^ 1|- Then (34) imphes rj^ ioi N > Nq via (36) and 
(37). We have" 

. N k-i \ / / ^ 

EE^(i+w) / E(i-w: 

^ fc=l 1=1 +^ / / \ ,.^1 

N k~l N-1 



r+(^) " (l-77/v)iV ^ 



(1 - w)iV ^ r+(*) (1 - 77A,)iV ^ r+(z) ' 
and thus (35) imphes iVei — >■ 0. 

To obtain Nej\j^i — > we can flip the process and consider Xj^ with the flipped rates as in the proof 
of Lemma 4.1. Now, in addition to (34), we also require the flip version of the first condition in (35): 

N-l - _ N-1 . N-1 

EN — i X - I\ — I X ^ I 

1=1 MO " h " it MO ^ ' 

which are guaranteed by the second condition in (35). □ 

Once we have the results of both Lemma 4.1 and Lemma 4.2, we can deduce that iVei — >■ and 
Ncn-i — > 0, which imply (32). Namely, from ei = cint^in + eioTTio, 

eiN < = ITT > 

since Lemma 4.1 ensures convergence of the denominator to 1 and Lemma 4.2 of the numerator to 
0. Similarly 

^ei iVei 

TVein < < s. 

°- ^10 - r_(0)/(r+(0)+r_(0)) 

since ttiq contains the positive probability (independent of A^) of an immediate return to 0. □ 

For a reaction system and splitting with unit net changes only, since splitting is unbiased we have 
Pi,i+i — Pi,i-i = ^, for i 7^ 0, N, and the contribution to r+(i) and r_(i) from splitting is N). 
Let us write ^{i,N) = ^{N)pi where ^{N) depends on N only (i.e., is state independent) and 
Pi — 0(1). Then, in any state, the contribution of the splitting is of 0(7(7V)), while the contribution 
of the reaction system is of 0{N) due to the standard scaling of reaction rates. Hence, we have the 
following result. 
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Theorem 4.1. If the reaction system has increments of size {1,-1} only; contains reactions aA — > 
{a-l)A+B, bB A+{b-l)B some a,b> 0; has rates with standard scaling K.f{N) = g'^bj^i-(a+b) . 
and if the splitting mechanism has increments of size {1, —1}, po,o = Pn,n = 1; with rate is j(i, N) = 
j(N)pi where j{N) and Pi = 0(1) satisfy: 

then the results of Proposition 4-2 apply with (3]^ = 1 and fpi = ^ -^-^^j^ ^ rio — X](a o -i)gi '^-'i ■ 
Proof. The transition rates for Xa are given by 

r+[i) = \-i{N)p, + N ~<\^lN)a.N{l-^IN\N. z = 0,...,iV-l, 

(a,6a)gl 

r_(z) = ^7(Ar)p, + iV ^ K!-\(i/iV),,jv(l-i/^)6,Ar, i = l,...,^. 

(a,&,-l)el 

On the boundary the rates are 

r+(0) = 7V r-_(iV) = 7V 

(o,b,i)ei (a,o,-i)ei 

and (33) holds with f+ = E(o,64)ei ' ^- = E(a,o,-i)ei '^-i- Also, 
- j{N)p, 

since k^'' > 0, where i? < oo is the number of reactions in the system. Therefore 

y |£jv(«)| < 2R max {^f }— r V - 
and the first condition in (38) ensures that J2i=i^ kA'(*)l ^ a^id (34) holds. On the other hand, 

Af-l . . N-1 . 

^ — A Z ^ — ^ 2 A ^ — A z 



and 



EN — I 2 x ^ JSI — I 



-[ r+{i) 7(7V) ^ 

so the last two conditions in (38) ensure that (35) are satisfied as well. □ 
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4.2 Example: Bistable Behaviour from Fast Splitting 



We revisit the same example of the reaction system we analyzed in Section 3.3: 

A B (39) 

b'^X A (40) 

A + B 2B (41) 

2A + b'X 3A (42) 

with the standard mass-action scaling, = A^^'^^+^^+^k"^. In this system the only reactions which 
counteract the absorption on the boundaries are the first two unimolecular reactions. Also, note 
that all system reactions change the molecular count of A only by increments of size 1. 

We chose the same simple splitting mechanism as before, since conditions (34), (35) or (38) are much 
easier to verify than conditions (31), (32). Recall that, if we were to assume j{N) — ^e^iV^ for 
some small > 0, then the limiting process for Xn would be the diffusion process in (28); the 
splitting noise is even less present if we were to assume 'f{N) = -^N, as shown in Section 3.3. In 
contrast, if we assume the rate 7(iV) grows fast enough so that A^^ In A^/7( A^) 0, then we can 
show that conditions in Proposition 4.2 are satisfied and the behaviour of the limiting process for 
Xn is described by a different two-state jump Markov process. 

S.01 „„j _ sio 



There are only two reactions in (39)- (42) active on the boundaries, so foi = and fio = To 

— fl - — ' 



verify (38) note that we have 7(«,iV) = -fiN)p, with = ^(1 - ^) = i{N - i)/N^, so 



^ ^ ~ ■^(N'\ ^ i(N - i\ ~ ^(N\ 2^ \ 7 



using partial fractions ^ ;^(t + {N-i) '^^^'^^ is the A^-th harmonic sum. Also 

7(A^) y K 7(iV) ^N-i 7(7V) 

and 



^ N ~i _ 1 1 



7(iV) ^ p. jiN) ^ z j{N) 
as well. Hence, A^^ln Af/7(A^) — > ensures that all conditions in (38) hold. 

This example shows that for any reaction system with unit increments whose drift has a double 
well potential, and for this particular choice of the splitting mechanism, we can identify orders of 
magnitude for 7(A^) that lead to different limiting behaviours: 

• If 7(Af) ^ A^, bistability is caused by large deviations of the Markov jump process, and the 
rescaled process transitions between neighbourhoods of the drift equilibirum points on a time-scale 
of order e^(T(^))"'*---2 , with N{-f{N))-'^ oo; 

• If j{N) ^ e^A^^, > a constant, bistability is caused by large deviations of a diffusion with 
a small perturbation coefficient, with transitions between neighbourhoods of the drift equilibirum 
points on a time-scale of order e'^ 

• If 7(A^) ^ A^'^ In A^, bistability is caused by excessive noise, and switching between the boundary 
points occurs on a time-scale of order 1. 
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Figure 3: Sample path Xjv(i) (left: x-axis=t, y-axis—XN{t) — N^^X{t) C [0,1]) of the system 
(22-25) with (27) splitting, under standard mass-action scaling for reactions and '^{N) = ^N^ 
(parameter N — 200); and the distribution of switching times plotted (dots) in terms of quantiles 
(right: x-axis=i, y-axis=fraction of switching times of length < t). Solid line (1 — e^*) indicates the 
quantiles of the exponential (mean 1) distribution for comparison. 

Note that the order of magnitude iV^ only represents the scale on which we have assumed that the 
variance of the splitting mechanism is in the diffusive case (see Assumption 7* in Sec 2.2). Also note 
that existence of two stable states in the deterministic model for the reaction system is not needed 
for the result of this Section. We chose the same reaction system in order to make the comparison 
with the results in Section 3 and emphasize the difference between the effects of "slow" and "fast" 
splitting on the same reaction system. 

Figure 3 shows an exact simulation of a sample path of the rescaled process Xj^ = N~^Xa for a 
relatively short period of time, spending most of its time at boundaries {0} U {!}, switching between 
them at approximately rates foi = '^i^ — l,^io = ^I'l = 1 (see Sec 3.3 for coefhcient values). 
Switching between states occurs at a time-scale Pfq = 1, and since fgi — — l,fio — — 1 
the distribution of switching times should approximately be an exponential distribution (mean 1) 
distribution. This is shown in the quantile plot in Figure 3, where the fraction of switching times of 
length < t is plotted against the same fraction 1 — e~* for the exponential (mean 1) distribution. 



5 Discussion 

We showed that there are two different types of stochastic bistable behaviour in which the system 
spends most of its time at or near one of two states and switches between them. For one of these 
types of bistability, because the magnitude of noise is high, it can occur even in a system whose 
deterministic model would not allow for a possibility of bistability at all. The detreministic system 
can have unique stable points, as for example in the neutral Wright-Fisher model with mutation. 
For the other type of bistability, where the noise is relatively low, one needs the reaction system to 
have two deterministic stable points, as for example in the Schlogl model. The important point is 
what constitutes "high" and "low" levels of noise: the determining quantity £a(A^) (H), depends 
on the relative size in terms of N of the variance to the average change in the system, where iV is a 
scaling parameter for the size of the system. We referred to £a{N) « as "slow" splitting, and to 
£a{N) w oo as "fast" splitting, interpreted relative to the reaction dynamics. 

We discussed the differences in the qualitative signatures of bistability in the two cases: 

• In case of "slow" splitting, the states where the process spends most of its time are determined by 
the drift of the deterministic model for the reaction system; in contrast, in case of "fast" splitting 



30 



they are simply the two extremes for the size ol the system. 

• In case of "slow" splitting, the rate of switching is determined by the relative magnitude of the 
splitting variance to the reaction drift and by the size of the potential barrier in the deterministic 
model for the system; on the other hand, in case of "fast" splitting, the rates of switching are 
determined only by the standardized rates of the reactions that are realizable from one of the 
extremes for the system size. 

• In case of "slow" splitting, the time-scale (3^ or /3j„ on which the switching happens is exponential 
in (some increasing function of) the size of the system; in contrast, in case of "fast" splitting, the 
time-scale jS^ is at most polynomial. 

We also showed that the observables of bistability (switching states and rates) are not sensitive 
to precise specification of the reaction system, as they depend only on: equilibrium points, size of 
potential barrier in "slow" splitting, and drift values at boundaries in "fast" splitting. However, 
bistability is very sensitive to the distributional form of the splitting/resampling mechanism: the 
variance of its distribution determines the potential barrier in "slow" splitting, and the harmonic 
sum of its transition probabilities determines the threshold for appearance of "fast" splitting. 

In the context of cellular systems of biochemical reactions, the problem of determining the partition- 
ing errors due to cell division is experimentally extremely challenging [22], [23]. The measurements 
for single cells rely on count estimates for related species rather than the molecular species of in- 
terest. In addition, in order to estimate the magnitude of intracellular noise, one has to separate 
the intrinsic from the extrinsic sources of randomness. How random is cell division, and how it 
compares in magnitude to the biochemical noise is a question that is very much open. However, 
since our analysis only depends on a few general features of the splitting mechanism (unbiasedness 
and time-scale of the rate), it is also possible that stochastic bistability is achieved by a set of auxil- 
iary reactions, instead of splitting, acting on a different time-scale from the rest of the system. For 
example, the protein bursting mechanism may act as the driver of stochastic bistability (Zong [_..], 
Kaufman [ ]). 

One can rely on the qualitative signatures of bistability in order to assess which of the two types of 
bistability we discussed is relevant in a specific cellular biochemical systems. For example, when the 
switching times are orders of magnitude greater than the molecular count of the switching species, as 
in the lysogenic switch of E. coli, the "slow" splitting is a more likely mechanism. In an experimental 
analysis of this system, [29] observed that the switching times of the cell are exponential in the 
number of protein burst events, and correspond to a calculation of the rare event probability of the 
bursts, as can be interpreted by large deviations in our "low" auxiliary noise ("slow" splitting) type 
of bistability. In contrast, when the switching times are relatively short, as in the gene expression 
switch in S. cerevisiae, the "fast" splitting is the probable mechanism. In the engineered chemical 
reaction network version of this system, [15] show that increasing the protein burst size (increasing 
the auxiliary noise) leads to more highly correlated switching behaviour in different cell lineages, as 
could be inferred from properties of our "high" auxiliary noise ("fast" splitting) bistability type. 
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